#include "element1d.h"
#include<cmath> //M_PI

//define PI M_PI

int main(int argc, char *argv[])
{
    const double PI = std::atan(1.0) * 4.0;
    const int N = 6;
    
    Point_1 mesh_pnts[N];
    mesh_pnts[0].set_coord(0.0);
    mesh_pnts[1].set_coord(PI / 4.0);
    mesh_pnts[2].set_coord(PI / 3.0);
    mesh_pnts[3].set_coord(PI / 2.0);
    mesh_pnts[4].set_coord(PI * 2.0 / 3.0);
    mesh_pnts[5].set_coord(PI);
    for (int i = 0; i < N; i++)
	mesh_pnts[i].set_index(i);
    for (int i = 0; i < N; i++)
	std::cout << mesh_pnts[i] << ", ";
    std::cout << "\b\b" << std::endl;
    Mesh_1 mesh(mesh_pnts, N);
    Quad quad;
    quad.load_data("interval.quad");
    quad.list();
    double A[N][N];
    double f[N];
    for (int i = 0; i < mesh.get_n_elements(); i++)
	std::cout << mesh.get_element(i);
    std::cout << std::endl;
    const Quad_Info& q_info = quad.get_quad_info(10);
    std::cout << "accuracy: " << q_info.get_accuracy() << std::endl;
    int n_pnts = q_info.get_n_q_points();
    std::cout << "no. of quad points: " << n_pnts << std::endl;
    std::cout << "points:" << std::endl;
    double cont = 0.0;
    for (int i = 0; i < n_pnts; i++)
    {
	const double *pAw = q_info.get_pAw(i);
	double x = quad.l2g(0, PI, pAw[0]);
	cont += x * x * x * x * x * x * x * x * x * pAw[1] * 2.0;
	std::cout << x << ", " << pAw[1] << std::endl;
    }
    double J = quad.Jacobian(0, PI);
    std::cout << "int: " << cont * J << std::endl; 
    return 0;
};

// Point_1 A;
// Point_1 B(0.3);
// Point_1 C(0.5, 1);

// Point_1 D;
// D = B;
// D.set_index(3);
    
// std::cout << B << D << std::endl;

// const double PI = std::atan(1.0) * 4.0; //Define a PI by myself, const means just one calculation.
//     std::cout << PI << std::endl;
    
//     const int N = 6;
//     Point_1 mesh_pnts[N];
//     //std::fstream mesh_file("mesh");
    
//     mesh_pnts[0].set_coord(0.0);
//     mesh_pnts[1].set_coord(PI/4.0);
//     mesh_pnts[2].set_coord(PI/3.0);
//     mesh_pnts[3].set_coord(PI/2.0);
//     mesh_pnts[4].set_coord(PI*2.0/3.0);
//     mesh_pnts[5].set_coord(PI);

//     for(int i = 0; i < N; i++)
//     {
// 	double pnt_coord;
// 	//mesh_file >> pnt_coord;
// 	mesh_pnts[i].set_index(i);
// 	//mesh_pnts[i].set_coord(pnt_coord);
//     }
//     for(int i = 0; i < N; i++)
// 	std::cout << mesh_pnts[i] <<std::endl;

//     std::cout << "\b\b" << std::endl;
//     //mesh_file.close();

//     Mesh_1 mesh(mesh_pnts, N);
//     int n_ele = mesh.get_n_elements();
//     for(int i = 0; i < n_ele; i++)
//     {
// 	Interval ele;
// 	ele = mesh.get_element(i);
// 	std::cout << ele << std::endl;
//     }
//     return 0;
